library(knitr)
# purl("h2fec_lifetime_MANOVA.Rmd", output = "h2fec_lifetime_MANOVA.R", documentation = 2)
rerun <- FALSE
set.seed(224819)
library(MCMCglmm)
## Loading required package: Matrix
## Loading required package: coda
## Loading required package: ape
library(tidyverse)
## ── Attaching packages ───────────────────────────────────────────────────────────── tidyverse 1.2.1 ──
## ✔ ggplot2 2.2.1 ✔ purrr 0.2.4
## ✔ tibble 1.4.2 ✔ dplyr 0.7.4
## ✔ tidyr 0.8.0 ✔ stringr 1.3.0
## ✔ readr 1.1.1 ✔ forcats 0.3.0
## ── Conflicts ──────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ tidyr::expand() masks Matrix::expand()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
library(parallel)
h2fec <- read.table("../../Data/Processed/eggs_per_vial.txt",
sep = '\t',
dec = ".", header = TRUE,
stringsAsFactors = FALSE)
# Standardize egg_total
h2fec$egg_total <- as.numeric(scale(h2fec$egg_total))
h2fec$animal <- seq(1, nrow(h2fec))
h2fec$treat <- as.factor(h2fec$treat)
pedigree <- h2fec[, c("animal", "sireid", "damid")]
names(pedigree) <- c("animal", "sire", "dam")
pedigree$animal <- as.character(pedigree$animal)
pedigree$sire <- as.character(pedigree$sire)
pedigree$dam <- as.character(pedigree$dam)
sires <- data.frame(animal = unique(pedigree$sire),
sire = NA, dam = NA, stringsAsFactors = FALSE)
dams <- data.frame(animal = unique(pedigree$dam),
sire = NA, dam = NA, stringsAsFactors = FALSE)
pedigree <- bind_rows(sires, dams, pedigree) %>%
as.data.frame()
genet_corr <- tibble(model = character(),
HS_STD = numeric(),
LY_STD = numeric(),
HS_LY = numeric(),
n_eff = numeric())
iter <- 6.5e6
burnin <- 5e4
thinning <- 500
chains <- 12
# 19 hours run time on nivalis (6e5 would be reasonable!)
if (rerun) {
HS <- h2fec %>%
filter(treat == "HS") %>% rename(Eggs_HS = egg_total) %>%
as.data.frame()
LY <- h2fec %>%
filter(treat == "LY") %>% rename(Eggs_LY = egg_total) %>%
as.data.frame()
STD <- h2fec %>%
filter(treat == "STD") %>% rename(Eggs_STD = egg_total) %>%
as.data.frame()
h2fec_mrg <- full_join(full_join(HS, LY), STD)
prior1 <- list(R = list(V = diag(3) * 1.002, nu = 1.002),
G = list(G1 = list(V = diag(3) * 1.002, nu = 0.002)))
priors <- list(prior1)
prior_names <- c("Tri: V = diag(3) * 1.002, nu = 0.02")
for (ii in 1:length(priors)) {
prior <- priors[[ii]]
fm <- mclapply(1:chains, function(i) {
MCMCglmm(cbind(Eggs_STD, Eggs_HS, Eggs_LY) ~ trait - 1,
random = ~ us(trait):animal,
rcov = ~ idh(trait):units,
data = h2fec_mrg,
prior = prior,
pedigree = pedigree,
family = rep("gaussian", 3),
nitt = iter,
burnin = burnin,
thin = thinning)
}, mc.cores = chains)
outfile <- paste0("../../Data/Processed/FEC_lifetime_model_prior", ii, ".Rda")
save(fm, file = outfile)
re <- lapply(fm, function(m) m$VCV)
re <- do.call(mcmc.list, re)
re <- as.mcmc(as.matrix(re))
n_eff <- effectiveSize(re)
# STD vs. HS
HS_STD <- re[ , "traitEggs_HS:traitEggs_STD.animal"] /
sqrt(re[ , "traitEggs_STD:traitEggs_STD.animal"] *
re[ , "traitEggs_HS:traitEggs_HS.animal"])
# STD vs. LY
LY_STD <- re[ , "traitEggs_LY:traitEggs_STD.animal"] /
sqrt(re[ , "traitEggs_STD:traitEggs_STD.animal"] *
re[ , "traitEggs_LY:traitEggs_LY.animal"])
# LY vs. HS
HS_LY <- re[ , "traitEggs_HS:traitEggs_LY.animal"] /
sqrt(re[ , "traitEggs_LY:traitEggs_LY.animal"] *
re[ , "traitEggs_HS:traitEggs_HS.animal"])
genet_corr <- bind_rows(genet_corr,
tibble(model = prior_names[[ii]],
HS_STD = median(HS_STD),
LY_STD = median(LY_STD),
HS_LY = median(HS_LY),
n_eff = mean(n_eff)))
}
write_csv(genet_corr, path = "../../Data/Processed/Genetic_Correlations_Fecundity.csv")
}
load("../../Data/Processed/FEC_lifetime_model_prior1.Rda")
fe <- lapply(fm, function(m) m$Sol)
fe <- do.call(mcmc.list, fe)
plot(fe, ask = FALSE)
plot(fe[, 1, drop = FALSE], ask = FALSE)
plot(fe[, 2, drop = FALSE], ask = FALSE)
plot(fe[, 3, drop = FALSE], ask = FALSE)
# Extract random effects, convert to mcmc.list
re <- lapply(fm, function(m) m$VCV)
re <- do.call(mcmc.list, re)
plot(re[, 1, drop = FALSE], ask = FALSE)
plot(re[, 2, drop = FALSE], ask = FALSE)
plot(re[, 3, drop = FALSE], ask = FALSE)
autocorr.diag(re)
## traitEggs_STD:traitEggs_STD.animal
## Lag 0 1.0000000000
## Lag 500 0.1832581246
## Lag 2500 0.0359859496
## Lag 5000 0.0073595947
## Lag 25000 0.0008738227
## traitEggs_HS:traitEggs_STD.animal
## Lag 0 1.000000000
## Lag 500 0.316749894
## Lag 2500 0.043083529
## Lag 5000 0.003399590
## Lag 25000 0.004914298
## traitEggs_LY:traitEggs_STD.animal
## Lag 0 1.000000000
## Lag 500 0.245468510
## Lag 2500 0.049304363
## Lag 5000 0.008819874
## Lag 25000 -0.001353121
## traitEggs_STD:traitEggs_HS.animal
## Lag 0 1.000000000
## Lag 500 0.316749894
## Lag 2500 0.043083529
## Lag 5000 0.003399590
## Lag 25000 0.004914298
## traitEggs_HS:traitEggs_HS.animal
## Lag 0 1.000000000
## Lag 500 0.293988216
## Lag 2500 0.035470069
## Lag 5000 0.005803579
## Lag 25000 -0.000953972
## traitEggs_LY:traitEggs_HS.animal
## Lag 0 1.000000000
## Lag 500 0.329957789
## Lag 2500 0.047373839
## Lag 5000 0.003177456
## Lag 25000 0.002173318
## traitEggs_STD:traitEggs_LY.animal
## Lag 0 1.000000000
## Lag 500 0.245468510
## Lag 2500 0.049304363
## Lag 5000 0.008819874
## Lag 25000 -0.001353121
## traitEggs_HS:traitEggs_LY.animal
## Lag 0 1.000000000
## Lag 500 0.329957789
## Lag 2500 0.047373839
## Lag 5000 0.003177456
## Lag 25000 0.002173318
## traitEggs_LY:traitEggs_LY.animal traitEggs_STD.units
## Lag 0 1.000000000 1.0000000000
## Lag 500 0.341618414 0.1646980549
## Lag 2500 0.070765607 0.0363956871
## Lag 5000 0.012885737 0.0049610233
## Lag 25000 -0.002303034 0.0006351579
## traitEggs_HS.units traitEggs_LY.units
## Lag 0 1.000000e+00 1.0000000000
## Lag 500 1.284304e-01 0.1628531527
## Lag 2500 1.603062e-02 0.0360801531
## Lag 5000 -1.557565e-05 0.0061217680
## Lag 25000 1.762078e-03 0.0002089429
effectiveSize(re)
## traitEggs_STD:traitEggs_STD.animal traitEggs_HS:traitEggs_STD.animal
## 73421.01 60571.97
## traitEggs_LY:traitEggs_STD.animal traitEggs_STD:traitEggs_HS.animal
## 64483.62 60571.97
## traitEggs_HS:traitEggs_HS.animal traitEggs_LY:traitEggs_HS.animal
## 66632.79 58873.16
## traitEggs_STD:traitEggs_LY.animal traitEggs_HS:traitEggs_LY.animal
## 64483.62 58873.16
## traitEggs_LY:traitEggs_LY.animal traitEggs_STD.units
## 52036.87 76920.44
## traitEggs_HS.units traitEggs_LY.units
## 98857.93 80116.07
# Concatenate samples
re <- as.mcmc(as.matrix(re))
head(re)
## Markov Chain Monte Carlo (MCMC) output:
## Start = 1
## End = 7
## Thinning interval = 1
## traitEggs_STD:traitEggs_STD.animal traitEggs_HS:traitEggs_STD.animal
## [1,] 0.24050995 0.061445782
## [2,] 0.25533472 0.023307804
## [3,] 0.10944533 0.023145763
## [4,] 0.26061087 0.075124197
## [5,] 0.32669058 0.130178556
## [6,] 0.17916471 0.073070329
## [7,] 0.09864812 0.002272214
## traitEggs_LY:traitEggs_STD.animal traitEggs_STD:traitEggs_HS.animal
## [1,] 0.08973243 0.061445782
## [2,] 0.16062404 0.023307804
## [3,] 0.04023973 0.023145763
## [4,] 0.06471364 0.075124197
## [5,] 0.21081884 0.130178556
## [6,] 0.24226351 0.073070329
## [7,] 0.10092543 0.002272214
## traitEggs_HS:traitEggs_HS.animal traitEggs_LY:traitEggs_HS.animal
## [1,] 0.03553152 0.01341190
## [2,] 0.02161892 0.02417345
## [3,] 0.04677608 0.04454500
## [4,] 0.06409154 0.04472212
## [5,] 0.05404178 0.08187758
## [6,] 0.03067786 0.09904113
## [7,] 0.02431111 0.01374877
## traitEggs_STD:traitEggs_LY.animal traitEggs_HS:traitEggs_LY.animal
## [1,] 0.08973243 0.01341190
## [2,] 0.16062404 0.02417345
## [3,] 0.04023973 0.04454500
## [4,] 0.06471364 0.04472212
## [5,] 0.21081884 0.08187758
## [6,] 0.24226351 0.09904113
## [7,] 0.10092543 0.01374877
## traitEggs_LY:traitEggs_LY.animal traitEggs_STD.units
## [1,] 0.04310317 0.4085918
## [2,] 0.10997286 0.3611172
## [3,] 0.04831385 0.6056401
## [4,] 0.05500526 0.6997679
## [5,] 0.16591677 0.3933644
## [6,] 0.33252455 0.4260540
## [7,] 0.11149806 0.7561495
## traitEggs_HS.units traitEggs_LY.units
## [1,] 0.2082286 0.4364787
## [2,] 0.2659669 0.4084736
## [3,] 0.2385426 0.3954161
## [4,] 0.1394668 0.4174038
## [5,] 0.2450056 0.3243062
## [6,] 0.2143837 0.2373222
## [7,] 0.2069289 0.3968859
# STD vs. HS
plot(re[ , "traitEggs_HS:traitEggs_STD.animal"])
plot(re[ , "traitEggs_STD:traitEggs_STD.animal"])
HS_STD <- re[ , "traitEggs_HS:traitEggs_STD.animal"] /
sqrt(re[ , "traitEggs_STD:traitEggs_STD.animal"] *
re[ , "traitEggs_HS:traitEggs_HS.animal"])
plot(HS_STD)
median(HS_STD)
## [1] 0.6838576
HPDinterval(HS_STD)
## lower upper
## var1 -0.5791482 0.9985265
## attr(,"Probability")
## [1] 0.95
# STD vs. LY
LY_STD <- re[ , "traitEggs_LY:traitEggs_STD.animal"] /
sqrt(re[ , "traitEggs_STD:traitEggs_STD.animal"] *
re[ , "traitEggs_LY:traitEggs_LY.animal"])
plot(LY_STD)
median(LY_STD)
## [1] 0.933697
HPDinterval(LY_STD)
## lower upper
## var1 0.3833676 0.999497
## attr(,"Probability")
## [1] 0.95
# LY vs. HS
HS_LY <- re[ , "traitEggs_HS:traitEggs_LY.animal"] /
sqrt(re[ , "traitEggs_LY:traitEggs_LY.animal"] *
re[ , "traitEggs_HS:traitEggs_HS.animal"])
plot(HS_LY)
median(HS_LY)
## [1] 0.7064098
HPDinterval(HS_LY)
## lower upper
## var1 -0.5952619 0.998148
## attr(,"Probability")
## [1] 0.95
library(tidyverse)
library(cowplot)
##
## Attaching package: 'cowplot'
## The following object is masked from 'package:ggplot2':
##
## ggsave
B <- data_frame(`HS vs. STD` = as.numeric(HS_STD),
`LY vs. STD` = as.numeric(LY_STD),
`HS vs. LY` = as.numeric(HS_LY))
colMeans(B)
## HS vs. STD LY vs. STD HS vs. LY
## 0.5101469 0.8425443 0.5199116
B %>% gather(Comparison, value) %>%
ggplot(aes(value, color = Comparison)) +
geom_line(stat = "density", size = 2) +
labs(x = "Genetic Correlation", y = "Density") +
theme(legend.position = c(0.15, 0.85),
text = element_text(size = 24),
legend.text = element_text(size = 18))
ggsave(last_plot(), file = "Genetic_Correlation_Plot.pdf",
width = 8, height = 6)
V_nu_to_a_b <- function(V, nu) {
require(tidyverse)
require(invgamma)
a <- nu / 2
b <- (nu * V) / 2
p <- tibble(
x = seq(0, 10, length = 200),
y = dinvgamma(x, a, b)) %>%
ggplot(aes(x, y)) +
geom_line()
print(p)
return(list(alpha = a, beta = b))
}
V_nu_to_a_b(1, 0.002)
## Loading required package: invgamma
## Warning: Removed 1 rows containing missing values (geom_path).
## $alpha
## [1] 0.001
##
## $beta
## [1] 0.001
if (rerun) {
prior1 <- list(R = list(V = 1, nu = 0.002),
G = list(G1 = list(V = diag(3) / 3, nu = 0.002)))
prior2 <- list(R = list(V = 1, nu = 0.002),
G = list(G1 = list(V = diag(3) * 0.5 * 0.7, nu = 0.002)))
M <- matrix(0.5 * 0.7, 3, 3)
diag(M) <- 0.5
prior3 <- list(R = list(V = 1, nu = 0.002),
G = list(G1 = list(V = M, nu = 0.002)))
priors <- list(prior1, prior2, prior3)
prior_names <- c("Sire: V = diag(3) / 3, nu = 0.002",
"Sire: V = diag(3) * 0.5 * 0.7, nu = 0.002",
"Sire: V = 0.5 diag, 0.35 off-diag, nu = 0.002")
h2fec$egg_total_std <- scale(h2fec$egg_total)
for (ii in 1:length(priors)) {
prior <- priors[[ii]]
fm <- mclapply(1:chains, function(i) {
MCMCglmm(egg_total_std ~ treat,
random = ~ us(treat):sireid,
data = h2fec,
prior = prior,
family = "gaussian",
nitt = iter,
burnin = burnin,
thin = thinning,
verbose = TRUE)
}, mc.cores = chains)
outfile <- paste0("sire_model_prior", ii, ".Rda")
save(fm, file = outfile)
re <- lapply(fm, function(m) m$VCV)
re <- do.call(mcmc.list, re)
re <- as.mcmc(as.matrix(re))
n_eff <- effectiveSize(re)
# STD vs. HS
HS_STD <- re[ , "treatHS:treatSTD.sireid"] /
sqrt(re[ , "treatHS:treatHS.sireid"] *
re[ , "treatSTD:treatSTD.sireid"])
median(HS_STD)
# STD vs. LY
LY_STD <- re[ , "treatSTD:treatLY.sireid"] /
sqrt(re[ , "treatSTD:treatSTD.sireid"] *
re[ , "treatLY:treatLY.sireid"])
median(LY_STD)
# LY vs. HS
HS_LY <- re[ , "treatHS:treatLY.sireid"] /
sqrt(re[ , "treatLY:treatLY.sireid"] *
re[ , "treatHS:treatHS.sireid"])
median(HS_LY)
genet_corr <- bind_rows(genet_corr,
tibble(model = prior_names[[ii]],
HS_STD = median(HS_STD),
LY_STD = median(LY_STD),
HS_LY = median(HS_LY),
n_eff = mean(n_eff)))
}
write_csv(genet_corr, path = "Genetic_Correlations.csv")
kable(genet_corr, digits = 3)
}
load("sire_model_prior1.Rda")
fe <- lapply(fm, function(m) m$Sol)
fe <- do.call(mcmc.list, fe)
plot(fe, ask = FALSE)
# Extract random effects, convert to mcmc.list
re <- lapply(fm, function(m) m$VCV)
re <- do.call(mcmc.list, re)
plot(re[, 1:3], ask = FALSE)
plot(re[, 4:6], ask = FALSE)
plot(re[, 7:9], ask = FALSE)
plot(re[, 10], ask = FALSE)
effectiveSize(re)
## treatHS:treatHS.sireid treatLY:treatHS.sireid treatSTD:treatHS.sireid
## 155811.5 153997.2 153740.5
## treatHS:treatLY.sireid treatLY:treatLY.sireid treatSTD:treatLY.sireid
## 153997.2 156405.7 157850.0
## treatHS:treatSTD.sireid treatLY:treatSTD.sireid treatSTD:treatSTD.sireid
## 153740.5 157850.0 155550.8
## units
## 154286.2
# Concatenate samples
re <- as.mcmc(as.matrix(re))
colnames(re)
## [1] "treatHS:treatHS.sireid" "treatLY:treatHS.sireid"
## [3] "treatSTD:treatHS.sireid" "treatHS:treatLY.sireid"
## [5] "treatLY:treatLY.sireid" "treatSTD:treatLY.sireid"
## [7] "treatHS:treatSTD.sireid" "treatLY:treatSTD.sireid"
## [9] "treatSTD:treatSTD.sireid" "units"
# ??? Heritability ????
HS <- re[, "treatHS:treatHS.sireid"] / (re[, "treatHS:treatHS.sireid"] + re[, "units"])
median(HS)
## [1] 0.02083006
LY <- re[, "treatLY:treatLY.sireid"] / (re[, "treatLY:treatLY.sireid"] + re[, "units"])
median(LY)
## [1] 0.0549402
STD <- re[, "treatSTD:treatSTD.sireid"] / (re[, "treatSTD:treatSTD.sireid"] + re[, "units"])
median(STD)
## [1] 0.3272483
# STD vs. HS
HS_STD <- re[ , "treatHS:treatSTD.sireid"] /
sqrt(re[ , "treatHS:treatHS.sireid"] *
re[ , "treatSTD:treatSTD.sireid"])
plot(HS_STD)
median(HS_STD)
## [1] 0.5223363
HPDinterval(HS_STD)
## lower upper
## var1 -0.9126065 0.9985669
## attr(,"Probability")
## [1] 0.95
# STD vs. LY
LY_STD <- re[ , "treatSTD:treatLY.sireid"] /
sqrt(re[ , "treatSTD:treatSTD.sireid"] *
re[ , "treatLY:treatLY.sireid"])
plot(LY_STD)
median(LY_STD)
## [1] 0.9045579
HPDinterval(LY_STD)
## lower upper
## var1 0.05970088 0.9995664
## attr(,"Probability")
## [1] 0.95
# LY vs. HS
HS_LY <- re[ , "treatHS:treatLY.sireid"] /
sqrt(re[ , "treatLY:treatLY.sireid"] *
re[ , "treatHS:treatHS.sireid"])
plot(HS_LY)
median(HS_LY)
## [1] 0.4888571
HPDinterval(HS_LY)
## lower upper
## var1 -0.8949345 0.9974266
## attr(,"Probability")
## [1] 0.95